## Measuring Political Support and Issue Ownership Using Endorsement Experiments,
## with Application to Militant Groups in Pakistan

library(R2jags)

odd <- function(x) x%%2==1

J <- 4 ## number of issue areas asked about
K <- 5 ## number of endorsement groups (including control)
n.province <- 4
n.division <- 26
province <- c(1,1,1,1,1,1,1,1,2,2,2,2,2,3,3,3,3,3,3,3,4,4,4,4,4,4)
n.per.division <- c(118,0,313,403,579,495,208,131,203,473,311,0,293,0,84,287,50,0,215,288,103,0,210,320,67,61) 
N <- sum(n.per.division) ## number of survey respondents
n.age <- 4

## Ideal Points
delta.division <- c(-1,-1,-0.75,-0.75,-0.5,-0.5,-0.25,-0.25, #province 1
		-0.5,-0.25,0,0.25,0.5,  #province 2
		-0.75,-0.5,-0.25,0,0.25,0.5,0.75,  #province 3
		0,0,0.5,0.5,1,1)  #province 4
division <- NULL
for (i in 1:n.division){  division <- c(division, rep(i,n.per.division[i])) }

lambda.division <- matrix(c(rep(0,n.division), rep(0, n.division), delta.division, -delta.division, delta.division/3), n.division, K)
omega <- c(0, 0.3, 0.3, 0.5, 0.5)

x <- rep(NA, N)
for (i in 1:N){ x[i] <- rnorm(1, mean=delta.division[division[i]], sd=1) }
x <- x / sd(x)

beta <- runif(J, 0.5, 2) ## discrimination parameter
alpha <- matrix(0.1,J,4)
for (j in 1:J){  ## variation in thresholds, but distance ensured
	while(alpha[j,2]-alpha[j,1] < 0.5 | alpha[j,3]-alpha[j,2] < 0.5 | alpha[j,4]-alpha[j,3] < 0.5){
		tmp <- rnorm(4, mean=0, sd=1)
		alpha[j,] <- tmp[order(tmp)]
	}}

s <- array(NA, c(N,J,K)) ##Individual level impact of group endorsement k
for (i in 1:N){
for (j in 1:J){
for (k in 1:K){
  s[i,j,k] <- rnorm(1, mean=lambda.division[division[i],k], sd=omega[k])
}}}

##Ordered Responses (5 choices)
Y <- array(NA, c(N,J,K))
for (i in 1:N){
  treat <- ifelse(odd(i),TRUE,FALSE) ## Half of Jake's respondents receive the control, remaining half split among treatments
for (j in 1:J){
for (k in 1:K){
  p1 <- pnorm(alpha[j,1] - beta[j]*(x[i] + s[i,j,k]))
  p2 <- pnorm(alpha[j,2] - beta[j]*(x[i] + s[i,j,k])) - (p1)
  p3 <- pnorm(alpha[j,3] - beta[j]*(x[i] + s[i,j,k])) - (p1 + p2)
  p4 <- pnorm(alpha[j,4] - beta[j]*(x[i] + s[i,j,k])) - (p1 + p2 + p3)
  p5 <- 1 - (p1 + p2 + p3 + p4)

  p <- c(p1,p2,p3,p4,p5)
  Y[i,j,k] <- sample(c(1:5), size=1, prob=p) ## Full array of potential outcomes
}}
if (treat==FALSE) { Y[i,,2:K] <- NA}
if (treat==TRUE) {
  Y[i,,1] <- NA 
  for (j in 1:J){ 
    answered <- sample(2:K,1); Y[i,j,-answered] <- NA  
}
}}

Y2 <- matrix(NA,N,J)
endorser <- matrix(NA, N,J)
for (i in 1:N){
	for (j in 1:J){
		endorser[i,j] <- (1:5)[!is.na(Y[i,j,])]
		Y2[i,j] <- Y[i,j,endorser[i,j]]
	}}
Y <- Y2


lambda.division.real <- lambda.division


#######################################
##  Prep and send data through JAGS  ##
#######################################
data <- list ("N", "J", "K", "Y", "endorser",
		"division", "n.division", "province", "n.province",
		"lambda.division.real") 
inits <- function() {list(alpha0=matrix(seq(-2,2,length.out=(4*J)),4,J), beta=runif(J, 0.2, 2), 
			x=rnorm(N), s=matrix(rep(0, N*J), N,J), 
			delta.division=rep(0,n.division), delta.province=rep(0,n.province),
			lambda.division=matrix(c(rep(0,n.division), rep(0,n.division*(K-1))), n.division, K), theta.province=matrix(c(rep(0,n.province), rep(0,n.province*(K-1))), n.province, K),
			omega=c(0,runif(K-1)) )} 
params <- c("lambda.division.error") 
model <- "ModelSimulationsDivision.txt" 
fit <- jags(data, inits, params, n.iter=20000, model.file=model)
fit <- autojags(fit, n.iter=20000, n.thin=20, Rhat=1.02, n.update=5) 
save.image(file = Outfile)


